#ifndef __NEWTONINTERP_H__
#define __NEWTONINTERP_H__

#include "HermiteInterp.h"

using namespace std;

// Newton interpolation 是 Hermite interpolation 的特殊情况，故可以继承
template <class T>
class NewtonInterp: public HermiteInterp<T,1>{
public:
  NewtonInterp();
  NewtonInterp(const InterpConditions<T,1>& IC);
  NewtonInterp(const string &_file);
  ~NewtonInterp();
  void addPoint(const T& x, const T& value);
  static T NevilleAitken(const T& x, const string &_file);
};

template <class T>
NewtonInterp<T>::NewtonInterp():HermiteInterp<T,1>() {};

template <class T>
NewtonInterp<T>::NewtonInterp(const InterpConditions<T,1>& IC):HermiteInterp<T,1>(IC){};

template <class T>
NewtonInterp<T>::NewtonInterp(const string &_file):HermiteInterp<T,1>(_file){};

template <class T>
NewtonInterp<T>::~NewtonInterp(){};

// 增加插值结点
template <class T>
void NewtonInterp<T>::addPoint(const T& x, const T& value){
  const int M = HermiteInterp<T,1>::tableOfDividedDiffs_.size();
  HermiteInterp<T,1>::tableOfDividedDiffs_[0].push_back(x);
  HermiteInterp<T,1>::tableOfDividedDiffs_[1].push_back(value);
  vector<T> vecn;
  for (int i = 2 ; i <= M ; i++){
       if ( i != M ){
            T tmp1 = HermiteInterp<T,1>::tableOfDividedDiffs_[i-1][M-i+1] - HermiteInterp<T,1>::tableOfDividedDiffs_[i-1][M-i];
            T tmp2 = HermiteInterp<T,1>::tableOfDividedDiffs_[0][M-1] - HermiteInterp<T,1>::tableOfDividedDiffs_[0][M-i];
            HermiteInterp<T,1>::tableOfDividedDiffs_[i].push_back(tmp1*1.0/tmp2);
       }
       else{
            T tmp1 = HermiteInterp<T,1>::tableOfDividedDiffs_[i-1][M-i+1] - HermiteInterp<T,1>::tableOfDividedDiffs_[i-1][M-i];
            T tmp2 = HermiteInterp<T,1>::tableOfDividedDiffs_[0][M-1] -HermiteInterp<T,1>::tableOfDividedDiffs_[0][M-i];
            HermiteInterp<T,1>::tableOfDividedDiffs_.push_back(vecn);
            HermiteInterp<T,1>::tableOfDividedDiffs_[i].push_back(tmp1*1.0/tmp2);
       }
  }
  vector<T> tmpvec;
  tmpvec.push_back(1);
  Polynomial<T> tmpPolynomial(tmpvec);
  Polynomial<T> vecPoly[M];
  vector<T> vec4Poly[M];
  for (int i = 0 ; i < M ; i++){
    if (i > 0){
      vec4Poly[i].push_back(-HermiteInterp<T,1>::tableOfDividedDiffs_[0][i - 1]);
      vec4Poly[i].push_back(1);
      vecPoly[i].setCoe(vec4Poly[i]);
      tmpPolynomial = (tmpPolynomial * vecPoly[i]);
    }
  }
  HermiteInterp<T,1>::interPoly_ = HermiteInterp<T,1>::interPoly_ + (HermiteInterp<T,1>::tableOfDividedDiffs_[M][0] * tmpPolynomial);
}

// Neville Aitken 加速算法
template <class T>
T NewtonInterp<T>::NevilleAitken(const T& x, const string &_file){
  InterpConditions<T,1> IC(_file);
  const vector<vector<T> >& table = IC.getTable();
  const vector<T>& vecx = table[0];
  const vector<T>& vecvalue = table[1];
  const int M = vecx.size();
  vector<T> vecdis;
  for (int i = 0 ; i < M ; i++)
    vecdis.push_back(x-vecx[i]);
  vector<T> ntable[M];
  for (int i = M - 1 ; i > 0 ; i--){
    ntable[i].resize(i,0);
    for (int j = 0 ; j < i ; j++){
         if (i == M - 1)
               ntable[i][j] = (vecdis[j] * vecvalue[j+1] - vecdis[j+M-i] * vecvalue[j])/(vecx[j+M-i] - vecx[j]);
         else
              ntable[i][j] = (vecdis[j] * ntable[i+1][j+1] - vecdis[j+M-i] * ntable[i+1][j])/(vecx[j+M-i] - vecx[j]);
     
    }
  }
  return ntable[1][0];
}


#else
//do nothing
#endif
